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Abstract 

We consider a dynamical model for a Fermi gas in the Bardeen-Cooper-Sclirieffer (BCS) superfluid state, trapped in 
a combination of a ID or 2D optical lattice (OL) and a tight parabolic potential acting in the transverse direction(s). 
The model is based on an equation for the order parameter (wave function) , which is derived from the energy density 
for the weakly coupled BCS superfluid. The equation includes a nonlinear self-repulsive term of power 7/3, which 
accounts for the Fermi pressure. Reducing the equation to the ID or 2D form, we construct families of stable ID 
and 2D gap solitons (GSs) by means of numerical simulations, which are guided by the variational approximation 
(VA). The GSs are, chiefly, compact objects trapped in a single cell of the OL potential. In the linear limit, the VA 
predicts almost exact positions of narrow Bloch bands that separate the semi-infinite and first finite gaps, as well as 
the first and second finite ones. Families of stable even and odd bound states of ID GSs are constructed too. We also 
demonstrate that the GS can be dragged without much distortion by an OL moving at a moderate velocity (~ 1 mm/s, 
in physical units). The predicted GSs contain ~ 10"^ — 10* and ~ 10^ atoms per ID and 2D settings, respectively. 
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1. Introduction 

Studies of matter-wave patterns in degenerate 
Bose [1] and Fermi [2] gases at ultra-low temper- 
atures have drawn a great deal of attention in the 
last decade. The theoretical description of patterns 
formed in a Bose-Einstein condensate (BEC) relics 
upon the Gross-Pitaevskii equation (GPE), which 
produces a remarkably accurate description of var- 
ious states, including solitons |3]. Experimentally, 
bright matter-wave solitons were created in BEC 
of ^Li [4] and ^^Rb atoms [5]. In either case, the 
interaction between the atoms was switched from 
repulsion to attraction by means of the magnetic 
field applied near the respective Feshbach resonance 
(FR) [6|7j. In BEC with repulsion between atoms, 
loaded in a periodic optical-lattice (OL) potential, 
which can be created as a standing wave by counter- 



propagating laser beams illuminating the conden- 
sate, localized states can be formed in the form of 
gap solitons (GSs). This was predicted taking into 
consideration the possibility of having a negative 
effective mass in parts of the bandgap spectrum 
induced by the periodic potential |8|9j . GSs are 
supported by the interplay of the negative effective 
mass and repulsive nonlinearity. They are stable 
against small perturbations [10| . The prediction was 
followed by the creation of a GS formed by ~ 250 
atoms of ®^Rb in a nearly one-dimensional (ID), 
i.e., "cigar-shaped", cross-beam optical trap, com- 
bined with the longitudinal OL potential |11I12) . 
In the experiment, the gas was subjected to accel- 
eration, in order to push atoms into a state with a 
negative effective mass. Other theoretically elabo- 
rated possibilities for the creation of GSs rely on the 
phase- imprinting method [13j . or on temporarily 
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adding a strong parabolic trap to the OL potential, 
with the objective to squeeze the condensate into a 
small region and then, relaxing the trap, to give the 
condensate a chance to self-trap into a compact GS 
state [14] . 

The creation of degenerate Fermi gases (DFG), 
and observation of the transition to Bardeen- 
Cooper-Schrieffer (BCS) superfluidity in them [15], 
suggest to consider possibilities of the creation of 
solitons in DFGs and BCS superfluids. Bright soli- 
tons were also predicted in Fermi-Bose mixtures 
with repulsion |16I17| or attraction [T7] between 
bosons and strong attraction between fermions and 
bosons. 

A possibility of the existence of fermionic GSs was 
mentioned in work [18] , in relation to a model of a 
binary boson condensate trapped in a one- or two- 
dimensional OL, with repulsion between the two 
components and zero intra-species interaction, as 
the respective system of coupled GPEs may also be 
realized as a system of equations for two fermionic 
wave functions, and thus the GSs reported in [18j 
may be interpreted as ID and 2D symbiotic two- 
component GSs in the binary DFG (term "symbi- 
otic" is frequently applied to two-component soli- 
tons in situations when each component in isola- 
tion is not able to form a soliton, while the coupling 
between them makes it possible [H]). Using more 
sophisticated models, one-dimensional GSs in the 
Fermi-Bose mixture with a small number of Bose 
atoms were predicted [20], and a possibility to cre- 
ate solitons in a mixed ID Bose-Fermi superfluid 
composed of a relatively small number of atoms was 
proposed in [5T]. 

In the strictly ID geometry, the DFG can be 
mapped into the Tanks- Girardeau gas of hard-core 
bosons [22], which, in a certain approximation, 
may be described by the ID nonlinear Schrodinger 
(NLS) equation with a repulsive quintic term [23] . 
In works [24] , GSs were predicted in the latter equa- 
tion which includes the OL potential, although the 
solitons reported in those works may be unstable 
against small perturbations. 

In this work, we aim to predict ID and 2D GSs 
in a BGS superfluid, which may be formed by Fermi 
atoms with weak attraction between ones with op- 
posite orientations of the atomic spin. The GSs will 
be found, chiefly, in the form of tightly bound states, 
localized in a single cell of the OL potential. A rig- 
orous many-body approach to the dynamics of the 
DFG, based on the system of quantum- mechanical 
equations of motion for individual fermions, made it 



possible to demonstrate the formation of fermionic 
GSs in a Bose-Fermi mixture [25] . However, such 
a microscopic approach becomes unfeasible as the 
number of fermions increases. A simplified dynam- 
ical model of the BCS superfluid relies on the use 
of a single wave function for fermions bound into 
Cooper pairs. This approach may be compared to 
that adopted in the Ginzburg-Landau theory, where 
the wave function describing the Fermi superfluid 
is introduced as a phenomenological complex order 
parameter [26]. Actually, the so developed descrip- 
tion of macroscopic objects, such as solitons, is valid 
in the hydrodynamic approximation, i.e., under the 
condition that the solitons' size is much larger than 
the de Broglie wavelength at the Fermi surface, in 
terms of the underlying microscopic distribution of 
the fermion atoms. As a result, the wave function 
of the BCS superfluid loaded in the ID or 2D OL, 
and confined in the remaining (transverse) direc- 
tion(s) by a tight "cigar-shaped" or "pancake" trap, 
obeys, respectively, the ID or 2D Schrodinger equa- 
tion, which includes the corresponding OL potential 
and a repulsive nonlinear term of power 7/3 [27] . 
The equation is derived from the Lagrangian which 
includes the energy density for the weakly-coupled 
BCS superfluid, as obtained in well-known works 
[28]. 

We construct GS solutions by means of numeri- 
cal computations, which are guided by an analytical 
variational approximation ( VA) |29| ; this approach 
has produced very accurate results in a model of the 
Bose-Fermi superfluid mixture [30] . The stability of 
the GSs is established in a numerical form, through 
direct simulations of the evolution of perturbed soli- 
tons. We conclude that the VA performs well in de- 
scribing profiles of compact fundamental GSs suffi- 
ciently deep in bandgaps, but it works poorly near 
bandgap edges, facing the difficulty in reproducing 
undulating tails demonstrated by numerical simula- 
tions in that case. In the linear limit, the approxi- 
mation very accurately predicts the location of the 
left edge of the first finite bandgap. 

The paper is organized as follows. The derivation 
of the ID and 2D equations for the wave function is 
presented in Section II. In Section III, we construct 
a family of stable fundamental GS solutions in the 
first two finite bandgaps of the ID model, using the 
VA based on the Gaussian ansatz and direct numer- 
ical solutions. The GS solutions maintain a tightly- 
bound (compact) shape, unless their chemical po- 
tential is taken very close to edges of the bandgaps. 
In Section III, we also present stable symmetric and 



antisymmetric bound states of fundamental GSs. 
and consider a possibility of dragging GSs by a mov- 
ing OL potential. In the framework of the GPE, the 
latter issue has drawn considerable attention as a 
means for transportation of matter-wave solitons, 
see work |31| and references therein. In Section IV, 
we report variational and numerical results for sta- 
ble GSs in the 2D equation. In addition to the case 
of the square-lattice potential, in Section IV we also 
find 2D radial GSs, supported by an axisymmet- 
ric potential which is periodic along the radius (in 
BEG trapped in the axisymmetric OL, radial GSs 
have been predicted in work [32]). Estimates for the 
number of atoms in the predicted ID and 2D GSs 
are given at the end of Sections III and IV, respec- 
tively, the result being A^iDsoiiton ^ 10"^ — 10'' and 

iV2Dsolito„=^ 2,000. 

For the sake of comparison, in Sections III and 
IV we additionally display variational and numerical 
results for GS families in the ordinary GPE-based 
ID and 2D models, which include the corresponding 
OL potential and, respectively, the cubic or quintic 
self-repulsive term. Actually, the accuracy provided 
by the VA in the model with the present model, with 
the nonlinearity of power 7/3, is higher than in the 
models with the cubic or quintic nonlinear terms. 

In addition to the spatially symmetric (even) fun- 
damental GSs, the ID equation gives rise to sub- 
fundamental (SF) solitons, which are localized odd 
states originating in the second bandgap. The SF 
solitons are also squeezed, essentially, into a single 
cell of the OL potential, but unlike the fundamen- 
tal GSs, they are (weakly) unstable. The SF solitons 
are considered in Appendix. In particular, the VA 
developed for them yields a physically relevant find- 
ing: in the linear limit, it accurately predicts the po- 
sition of the narrow Bloch band separating the first 
and second bandgaps. 

2. Dynamical Equations 

2.1. The equation in the general form 

We consider a BCS superfluid of Fermi atoms 
with spin 1/2 and mass m, assuming weak attraction 
between fermions with opposite orientations of the 
spin. To derive an effective equation for the super- 
fluid order parameter, we start with the well-known 
expression for the energy density of the weakly cou- 
pled BCS superfluid in the 3D space |28I33| , 

SD = (3/5)p3D£f + 4TTah^plj^/i2m) + ..., (1) 



which is built as an expansion in powers of kpa, 
where hkp and a are the Fermi momentum and 
scattering length of the weak interaction between 
fermions in the BCS limit, ep — (hkp) /{2m) is 
the Fermi energy, and p^D the atomic density. The 
second term in Eq. ([1]) is a small correction to the 
first one due to the underlying condition, kpa <ti 1 
(which actually implies that a is much smaller than 
the de Broglie wavelength at the Fermi surface). 
Taking into regard the Cooper pairing of spin-up 
and spin-down fermions, their total density is p^^ = 

2{27t)-^ J^^ Airk'^dk = {37T^)~' {2meF/h^Y^\ 
From here, the Fermi energy can be expressed in 
terms of the atom density. 



^^ fo 2 \2/3 



(2) 



hence the energy density in Eq. ([T]) is cast in the 
form of 



„ _ 3{3n^)y'h^ 5/3 ^ 2nah^ ^ 



(3) 



which includes the lowest-order correction propor- 
tional to scattering length a. 

In the framework of the density- functional theory, 
the BCS superfluid is described by a complex order 
parameter (wave function) , ^, such that p^j~, = |^p. 
If the Fermi energy is much larger than the depth of 
the external potential, V{r), the evolution of the or- 
der parameter can be derived from the correspond- 
ing Lagrangian density which includes energy den- 
sity dS]), 



£ = —(***( -***) 



|VvI'|2_y(r)|M'|2 

2moff 

3(3^^)^l^|io/3_2rf|^|4^ 

10?n m 



(4) 



where rUcs is the effective mass of the order- 
parameter field in the density-functional theory, 
which determines the gradient term in the La- 
grangian density. Such a formalism has been used 
in the description of Fermi superfluids in various 
contexts [51] . 

The Euler-Lagrange equation which follows from 
the full Lagrangian, L = J Cdr, where density ([4]) 
is inserted, takes the form of the three-dimensional 
NLS equation with the repulsive nonlinear term of 
power 7/3, and a small correction in the form of the 
cubic term: 



£: 



ih-^t 



2moff 



(9!: 



dl + af) * 



2m 



>'f"m 



4/3 



87ra|*r 



1' + y(r)*. (5) 



Since \'$\'^ is defined as the atomic density, Eq. (O 
is supplemented by the normalization condition, 



jjjmrMd.dyd. 



^N, 



(6) 



where N is the total number of atoms. 

Equation ([S]), which was derived from the local 
Fermi distribution, applies to the description of 
spatially nonuniform patterns in the hydrodynamic 
limit, which assumes that the nonuniformity docs 
not strongly disturb the local distribution. To say 
it more accurately, the hydrodynamic approach is 
valid if the characteristic size of the macroscopic 
pattern (actually, the width of the soliton, which 
is close to the OL period, A/2, sec Fig. [2] below) is 
much larger than the de Broglie wavelength at the 
Fermi surface, i.e.. 



A > 27r/fc/ 



(7) 



This condition is similar to that necessary for the 
validity of the hydrodynamic approximation in clas- 
sical rarefied gases: the scale of the flow must be 
much larger that the free-path length. 

2.2. The one- dimensional equation 

Taking V{r) as the potential of the 3D OL, Eq. 
([5]) can be used to construct three-dimensional 
GSs. However, producing systematic results for 3D 
solitons is a challenging numerical problem. In this 
work, our objective is to reduce the equation to its 
ID and 2D forms, assuming, as mentioned above, 
that the superfluid is held in a relatively tight cigar- 
shaped or pancake-like trap, that correspond to the 
following transverse potentials: 



Vi^'^^ = imc^i (y' 



z^) ■ yfp""^) 



1 2 2 

—muj I z . 
2 



(8) 

For this reduction to be valid, the Fermi energy must 
be larger than both the distance between energy lev- 
els in the spectra of harmonic-oscillator potentials 
(III) , and depth Vq of the longitudinal OL potential: 

ep > fiuj±,VQ . (9) 

Equations (O and ([9]) constitute conditions neces- 
sary for the validity of the approach elaborated in 
the present work. 

The next step is reducing Eq. ([5]) to effective ID 
and 2D equations. In the ordinary GPE with the cu- 
bic nonlinearity, the reduction of the 3D equation to 



its ID and 2D forms can be performed in different 
ways, depending on the particular setting |35I36I37] . 
In the simplest situation, the reduction of the 3D 
equation to ID starts with assuming the factoriza- 
tion of the wave function. 



^(x, y, z, t) = ^{x, t) exp 



y'^ + z"^ 



(10) 



where the transverse harmonic-oscillator length is 
Oho = \Jf^/ {'m^±)i and, in the case of a very tight 
confinement, n = 1. However, due to condition ([5]), 
£p corresponds not to the ground state of the trans- 
verse potential, but rather to an excited state with 
large quantum number n, so that ep = nhLu±, and 
the transverse size of the domain filled by the Fermi 
superfluid is ~ v^aho- Accordingly, the second mul- 
tiplier in Eq. (|10[) approximates a lumped superpo- 
sition of excited states of the transverse harmonic 
oscillator, with the corresponding quantum number 
taking values from 1 to n. The self-consistency de- 
mands that Ep = nhLu± , which, after a simple anal- 
ysis, leads to relation 

2n = {{3/2)7:^ p,^aLf'. (11) 

One-dimensional function ^(x,t) in ansatz ()10p 
accounts for the dynamics in the x direction, its nor- 
malization being determined by Eq. ([6]). The substi- 
tution of Eq. (fTO|) in Eq. ([5]) and subsequent averag- 
ing of the 3D equation in the transverse plane yield 
the effective ID equation 



ih^t^-- 



2moff 
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X $, 



(12) 



where the ID potential with strength — e corre- 
sponds to the OL with period A/2 (the negative 
sign in front of e implies that a local minimum of 
the potential is set at a; = 0, where the center of the 
soliton will be placed). 
By means of rescaling 



2N , , TOcffA ~ A , 

* = V ^-"ho V, t = t, X = —X, 

V nX Att^Ti 27r 

Eq. p^ is cast in the dimcnsionlcss form. 



(13) 



+.9id| 



^V' " Vo cos (2a;) ip 



(14) 



(tildes are dropped here), with the rescaled ID wave 
function subject to normalization 



+ CXD 



\^{x)\Ux==l. 



(15) 



The effective strengths of the Fermi nonhncarity, 
weak cubic interaction, and OL potential are defined 
here as fohows: 



G 



(7/3) ^ 3 mcff / 3A^7V ^ 



ID 



10 m \8Tma^^ 

2 



Vb = mc 



A 
2Trh 



e, 



TOcff aAiV 
5iD = ^ o_„,2 (16) 



771 27r77a 



ho 



ActuaUy, Vq is the OL depth measured in units of 
the recoil energy, Eji — {2nhy / {2mX^) . Below, it 



will be demonstrated that G 



(7/3) 
ID 



takes values < 10 



(see Fig. [5]), while gm is confined to the range of ~ 
10"'^, hence the cubic term may be safely neglected 
in the first approximation. 

2.3. The two-dimensional equation 

Under the transverse confinement in direction z, 
which corresponds to the "pancake" configuration, 
Eq. ([5]) can be reduced to a 2D form in the plane 
of {x,y). To this end, following [3S], we substitute 
^{x,y,z,t) = $(x,?/,i)exp(-zV(277a^J), cf. Eq. 
(jlOp . Averaging in z and making use of the same 
rescalings as in Eq. p3p , except for a different trans- 
formation of the wave function, 

2^3/4^ 



77 ' Ay/aiio 

we arrive at the following 2D equation (tildes are 
again dropped here) : 

+32d|^|'V' - Vo [cos (2x) + cos (2y)] 7/>, (18) 

where we have included the 2D lattice potential and 
adopted the following definitions: 
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m flho 
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1. 



(19) 



(20) 



Finally, the above-mentioned self-consistency con- 
dition, nhuj± ~ ep, leads, in the present case, to 
relation 



2?7 



V2 



'^'^Pwaio 



2/3 



(21) 



cf. its counterpart (fTTj) for the cigar-shaped config- 
uration. 

2.4. Other equations 

The derivation presented above should be modi- 
fied in the case of extremely tight transverse confine- 
ment, with ep ^ huj^. If the very tight trap is cigar- 
shaped, the derivation must start from the ID Fermi 
distribution, with the respective atomic density re- 
lated to ep as follows: pijy = 2(27r) /_fc^ dk = 
{2/TTh) y/2mep, hence 

ep = ir^h^plo/iSm). (22) 

It is known that the energy density of the ID Fermi 
superfluid is [38] £id = (1/3)pid£f- Using expres- 
sion (im, one can derive the energy density of the 
ID superfluid [39]: 

iu^{7:^fiy24m)pl^. (23) 



£ 



Similarly, for the 2D superfluid subjected to an ex- 
tremely tight confinement in direction z, the energy 
density is £20 = {^/'^)P2D^f, the 2D atomic density 
being p2D = 2(27r)^^ /^^ ^ 27rfcdfc = [m/nh^) ep, 
i.e., ep = Trh^p2£,/m. Thus, the energy density of 
the 2D superfluid can be written in terms of the 2D 
density [lU] . 

S2D = (7rftV2777) pI^ . (24) 

Repeating the analysis which produced Eq. ([5]), 
we arrive at the following 2D and ID equations cor- 
responding to the limit case of the very tight trans- 
verse trap: 






+V{r)'iP 



(25) 



The equation with the cubic nonlinearity, in the first 
line, is formally equivalent to the ordinary GPE in 
two dimensions, while the ID equation in the sec- 
ond line coincides with the above-mentioned quintic 
equation for the ID gas of hard-core bosons [23]. 

Assuming that V represents the OL potential 
in the ID variant of Eq. ([25]) . the respective fam- 
ilies of GS solutions [27] are mathematically tan- 
tamount to the ID bosonic GSs generated by the 
above-mentioned NLS equation with the quintic 
self- repulsive term [53]. However, unlike Eq. (fH|) . 
the physical relevance of the ID variant of Eq. (|25)) 
for the prediction of GSs in the BCS superfluid is 



impugnable, because subsequent analysis demon- 
strate that such solitons, on the contrary to those 
produced by Eq. p4|) . contain few (^ 100) atoms 
[27] , which makes the very concept of the superfluid 
doubtful in such a situation. Similarly, it is possible 
to construct a family of GSs in the 2D version of Eq. 
(P5)l . assuming that V{r) is the 2D OL potential. 
Mathematically, they will be equivalent to the GS 
solutions of the ordinary 2D GPE [8]. However, in 
this case too, the number of atoms in such 2D soli- 
tons, if they are realized in terms of the BCS super- 
fluid, turns out to be small (on the contrary to the 
2D solitons predicted by Eq. ([18])). Therefore, we 
will focus on the physically relevant models based 
on Eqs. (fT4|l and (|18|) . The results will be compared 
to those for the bosonic GSs generated by one- and 
two-dimensional GPEs, (1551) and (H51). see below. 



3. One-dimensional solitons 

3.1. Variational approximation 

Stationary solutions to Eq. (fT4|) are looked for 
in the usual form, ^jj{x,t) = e~'^*0(j;), where /i is 
the chemical potential, and real function (/)(a;) obeys 
equation 

^ic^+{l/2)f-G''^'^W^/'^ + Vocosi2x)(l) = 0, (26) 

with (f) = (P4'/dx'^, in which the small cubic term 
is dropped (its effect will be considered below) . This 
equation, together with normalization condition 
rS)) . can be derived from the following Lagrangian, 



L 



+ OC 



/^r 



:G 



(7/3) ,10/3 
ID V 



+Vq cos{2x)(I)'^] dx — iJ., 



(27) 



by demanding 5L/54) = dL/df-L = 0. To apply the 
VA, we use the Gaussian ansatz l29l. 



<j>{x) 



-1/4 M 



2W^ 



(28) 



where variational parameters are the soliton's norm 
and width, M and W (in addition to /i). The use 
of the Gaussian is justified by the numerical re- 
sults showing that, except for a narrow vicinity of 
bandgap edges, the fundamental GSs are strongly lo- 
calized solutions, see below. We also note that ansatz 
(PS)) implies that the center of the soliton is placed at 
a local minimum of the OL potential. A straightfor- 
ward generalization of the VA, which treats the cen- 
tral coordinate as another degree of freedom of the 



GS, demonstrates that solitons may also be found 
with the center located at a local potential maxi- 
mum, but they are always unstable against the shift 
from that position. 

The substitution of ansatz (|28p in Lagrangian ([27]) 
yields 



i = /i(A/-l)- 



M 
4VK2 



VoMe 



-w^ 



3^3/2^(7A3)^^ 

7rV6 1^2/3 • 



(29) 



The first variational equation following from Eq. 
([29)) . dL/dfi = 0, yields (as it should) M = 1, there- 
fore, M = 1 is substituted in other variational equa- 
tions below, except for equation dL/dM = 0, where 
M = 1 is substituted after the differentiation. The 
remaining variational equations are dL/dW = 0, 
which predicts a relation between the soliton's width 



and effective nonlinearity strength, G 



(7/3) 
ID 



4V3G 



53/2,^1/3 



(7/3) 



(30) 



and dL/dM = 0, which yields ^ as a function of W 
and G 



(7/3). 
ID 



fJ- 



vsGr^'^ 



"^ID 
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Voe 
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(31) 



The first noteworthy consequence of the above 
equations is that, for Vb not too small, they predict a 
certain value of ^ at G\-^ — > 0, which corresponds 
to the stationary linear Schrodinger equation. In- 
deed, setting gJd^^^ = in Eq. ^ yields 



4Vb = W^~''exp(W^2^. 
This expression attains a minimum, 

(^o)„,i„ = eVl6 « 0.462, 



(32) 



(33) 



at W^ = 2. At Va > (Vo)nii„, Eq. ([32]) yields two 
solutions for W^ . 

A commonly known fact is that the linear 
Schrodinger equation with the periodic potential 
cannot give rise to localized states. Nevertheless, 
the results presented in Fig. [U which displays the 
well-known band structure for the linear version of 
Eq. ([26]) . together with the /i (Vb) curve (the con- 
tinuous red line), as obtained from Eqs. ([32]) and 
([5T]) with G\j^ = 0, clearly indicate that, if Vq ex- 
ceeds minimum value p3p . the variational solution 
makes sense in the linear limit: it accurately pre- 
dicts the location of the left/lower edge of the first 




Fig. 1. (Color online) The bandgap structure of the linearized 
version of Eq. I|26|l . Narrow regions between the gaps repre- 
sent Bloch bands. The solid and dotted curves are borders 
between the semi-infinite and first finite gaps, and between 
the first and second finite gaps, as predicted by the varia- 
tional approximation (see text). 

finite bandgap, as tlic starting point of the GS fam- 
ily inside this bandgap, see Fig. [5] below. In other 
words. Fig. [1] shows that the VA makes it possible 
to predict the location of the narrow Bloch band 
which separates the semi-infinite gap and first finite 
bandgap (a similar result was obtained by means of 
VA in [H]). In particular, in the case of Vb = 5, for 
which most results are presented below, Eqs. ([5^ 
and (PT|) with G\-^ = predict the left edge of the 
first finite bandgap at 



(.. ^(var) 
1^1 /left 



("^^0 = 5) 



(34) 



This value almost exactly coincides with its coun- 
terpart found from the numerical solution of Math- 
ieu equation (i.e., the linearization of Eq. (|26[) . 
K)L"r^ (^0 = 5) « -2.893. ^ 

The above-mentioned variational curve shown in 
Fig. [T] corresponds to a smaller root of Eq. ([5^ . 
Through the other (larger) root of the same equa- 
tion, the VA is actually trying to predict another 
narrow Bloch band, which separates the first and 
second finite bandgaps in Fig. [TJ This root produces 
a large error in predicting the border between the 
first and second gaps, as the underlying even ansatz, 
adopted in Eq. (pS)) . is inadequate in that case. How- 
ever, the border can be accurately predicted by a 
modified version of the VA, based on a properly 
modified (odd) ansatz, see Eq. (|45p in Appendix and 
the dotted blue curve in Fig. [1] 

In addition to the VA based on Gaussian 
ansatz ([25)) . we have also elaborated its coun- 
terpart based on the hyperbolic secant, (l){x) = 
y/M/2Wscch{x/W). Eventual results produced 



by this ansatz, which we do not display here, are 
quite close to those generated by the Gaussian, al- 
though the accuracy is slightly lower; in particular, 

it predicts (/^i)LT^ (^o = 5) - 2.854, cf. value ^ 
produced by the Gaussian-based VA. 



3.2. Numerical results: fundamental gap solitons 

Numerical solutions where obtained by inte- 
gration of time-dependent Eq. (jl4|) (and time- 
dependent counterparts of bosonic Eqs. ([35]) . see 
below), using the Crank-Nicholson discretization 
scheme, until the solution would converge to a 
time-independent real form. The equations were 
discretized using time step 0.0005 and space step 
0.025, in domain -20 < a; < +20. This method of 
obtaining the stationary solutions guarantees that 
they are stable. The results presented here were 
obtained dropping the small cubic term in Eq. (|14p : 
effects induced by this additional term are consid- 
ered below separately. 

The family of fundamental GSs predicted by the 
VA is characterized by dependence G\j^ (//), which 
was obtained from a numerical solution of Eqs. (|30p 
and (j31[) (it can be easily translated into a depen- 
dence between fi and the number of atoms, using 
Eqs. (fT6|) and (fTTjl ). This dependence is displayed for 
two different values of the OL strength, Vb, in Fig. 
[51 For the sake of comparison, the figure includes 
similar dependences for families of ID bosonic GSs, 
as obtained by means of the VA, based on the same 
ansatz ([55| . from the following bosonic equations 
with the cubic and quintic self-repulsive terms. 



M0 + -^(f)" - 



(3) -3 
(5) ,5 



+ Vq cos(2x)(t) = 0. (35) 



Here, wave function (t>{x) is subject to normalization 
condition ^, and G^^ = (Xajnal^) N [5155] . 



G 



(5) 
ID 



72) N^ 



with N the number of bo- 



son atoms, tts the scattering length characterizing 
the repulsive interactions between them, and aho the 
same transverse-trapping size as above. 

Comparison of typical shapes of the numerically 
found GSs with the respective VA profiles predicted 
by the Gaussian ansatz is displayed in Fig.[3|(a set of 
GS shapes obtained from Eq. ([35]) with the quintic 
nonlinearity, which is also displayed in this figure, 
makes it possible to compare the GSs in the BCS su- 
perfluid with their bosonic counterparts). The figure 
includes examples of the GSs belonging to both the 
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Fig. 2. (Color online) The nonlincarity coefficient, G^^ 
Eq. l|14|l for the BCS superfluid, versus chemical potential /^, 
for the family of fundamental ID gap solitons found in the 
first two finite bandgaps of periodic potential — Vbcos(2x), 
with Vb = 5 (a) and Vb = 3 (b). For comparison, also shown 
are dependences between the respective nonlinear coefficient, 



G*?' or G 



(5) 



and /i for families of fundamental gap solitons 



in two ID bosonic equations 1351 1. Here and in similar figures 
displayed below, shaded vertical areas represent Bloch bands 
which separate the gaps. Numerically found solution families 
are depicted by continuous curves, whereas the predictions 
produced by the variational approximation are shown by 
chains of symbols. 



first and second bandgaps, which can be easily iden- 
tified by values of the chemical potential indicated 
in the panels. In addition to the numerical and vari- 
ational profiles, in Fig. [3] we also plot their simplest 
counterparts predicted by the Thomas-Fermi (TF) 
approximation, which were obtained, as usual [1], 
by dropping the second-derivative term in Eqs. (|26p 
or (|35[) (the TF profiles are not shown for very weak 
nonlinearity. where this approximation is irrelevant) 
Figure[3]demonstrates that the fundamental GSs. 
unless taken near edges of the bandgaps, are in- 
deed compact objects with a quasi- Gaussian shape, 
trapped in a single cell of the OL potential. This 
shape radically changes as one comes very close to 
a bandgap's edge, or pushes the solutions to higher 
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Fig. 3. (Color online) One-dimensional GS wave forms (nor- 
malized to unity) in the first and second bandgaps as ob- 
tained (a) from Eq. 1261 1. without the small cubic term, and 
(b) from bosonic Eq. 1351 1 with the quintic nonlinearity, for 
Vb = 5. Here and below, labels "num", "var", and "TF" 
represent numerical, variational, and Thomas-Fermi results, 
respectively and the thin sinusoidal line depicts the under- 
lying periodic potential, — Vbcos(2a::). 



values of /i (in particular, the curves representing 
the GS family in Eq. (j35p with the cubic nonlinear- 
ity are aborted in Fig.[2](b) at a spot where the soli- 
ton undergoes a transition to a complex shape with 
undulating tails, making the further use of the VA 
irrelevant). 

The variational solutions displayed in Fig. [UJa) 
emerge, at G = 0, at the value of /.t given by Eq. 
([M|) (recall that it almost exactly coincides with the 
actual left edge of the first bandgap; the same is 
true for Fig.I^l^b). However, an essential defect of the 
variational solutions for these families is that they 
ignore the Bloch band separating the first and sec- 
ond bandgaps, going across it (inside the bands the 
VA predicts the so-called quasi-solitons, i.e., nearly 
localized solutions, with the smallest amplitude of 
nonvanishing tails attached to the central "body" 
|42|). Another noteworthy inference is that, as seen 
from Figs. [2] and [31 the accuracy provided by the 



VA for the description of the GS family generated 
by Eq. (p6|) for the BCS superfluid is better than for 
the bosonic GS families. 

3.3. Physical estimates for the fundamental gap 
solitons 

Undoing rescalings (fT5|) and P^ . we conclude 
that conditions (O and ^. which underlie the 
derivation of the effective ID Eq. p^ . are definitely 
satisfied for G\^ > 3, i.e., according to Fig. [2l 
starting from the middle of the first finite bandgap. 
To assess the feasibility of the creation of the GSs 
in the BCS superfluid trapped in the OL, it is nec- 
essary to estimate the expected number of atoms, 
N, in the predicted quasi-lD soliton. Getting back 
to the physical units, we arrive at the following 
estimates for N and the corresponding value of n 
(recall it is the largest quantum number of the filled 
states in the transverse trapping potential, which 
appears in Eqs. pO)) and plj) ): 



N- 



IDsoliton ■ 



10^- 



"Icff 



3/2 



G 



■5/2 



T) («ho/A)^ 



n^ ~ 10 (iViDsoiitonaho/A)^ 



(36) 



Figure [5] demonstrates that the largest achiev- 
able value of the effective nonlinearity strength is 

(^iD ) ^10. Then, assuming Oho/A ~ 1 (for 

V / max 

^Li atoms, this implies the use of trapping frequency 
uj± ~ 27r X 1 KHz, if the OL is induced by light 
with wavelength A ~ 1 //m), relations p6p predict 
-^IDsoliton in the range of 10^, with n ~ 100. For a 
larger wavelength, A ~ 3 /im, one concludes that 
-^IDsoliton takes values in the range of 10^, with n 
~20. 

3.4. Bound states of gap solitons 

We have also constructed symmetric and antisym- 
metric ("twisted") bound states of the fundamen- 
tal GSs (unlike the SF solutions, considered in Ap- 
pendix, they are stable solutions to Eq. (|14p ). These 
states were formed by placing in-phase or out-of- 
phase pairs of identical GSs (ones with equal or op- 
posite signs, respectively) in adjacent cells of the 
OL potential and allowing them to achieve a stable 
configuration. Typical examples are shown in Fig. |3] 
(a). In the experiment, the relative phase of adjacent 
atomic clusters may be controlled by laser beams. 




Fig. 4. (Color online) (a) Examples of symmetric and anti- 
symmetric ("twisted") stable bound states of two ID funda- 
mental GSs, obtained from Eq. I I14I I. (b) Numerically con- 
structed (continuous lines) plots of the nonlinearity versus 
the chemical potential for symmetric (S) and antisymmetric 
(A) bound states, for the BCS superfluid and two bosonic 
models. The analytical results (chains of symbols) are gen- 
crated using Fig. [J] with Gid rcscaled as described in the 
text. 



The Gid(m) curves for the synmretric and an- 
tisymmetric bound states are shown in Fig. [41(b) 
(again, they are presented together with similar de- 
pendences for bound GS pairs generated by two 
bosonic equations ([55)1 ). As might be expected, these 
curves can be obtained, in an almost exact form, 
from those for the fundamental solitons (see Fig. 
El), by means of rescalings: G^d^^^ -^ 2^^^G?^^^ 



(3) 
ID 



(3) ^(5) 






.(5) 



'ID 



ID 



4:G\j^, respectively. Indeed, 



G 

if the equations are written in the notation with 
a fixed nonlinearity coefficient and variable norm, 
the rescalings simply imply that the norm of the 
bound state is twice that of the fundamental soliton. 
The fact that the so defined norm of the symmetric 
bound states slightly exceeds the norm of their an- 
tisymmetric counterparts is natural too, as the van- 
ishing of the density at the center of the antisym- 
metric state makes its total norm slightly smaller. 
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represented in Fig. [5] The above physical estimates 
suggest that typical length and time units in the 
present setting are ~ 1 /Ltm and 1 ms, from which 
we infer that the stable transportation of the GSs is 
possible for velocities < 1 mm/s. 
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Table 1. The relative loss of atoms suffered by the 
one-dimensional gap solitons at the end of dragging 
shown in Fig. [5l 

4. Two-dimensional solitons 

4.1. Variational approximation 
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Fig. 5. (Color online) Dragging GSs with G-!„ — ■^-.^ 
by the OL potential — 5cos[2(a; + vt)\, which was suddenly 
set in motion at constant velocities —v = —2, —4, —8. The 
upper (BEG) and lower (BGS superfluid) panels display a set 
of snapshots of the density, taken at consecutive moments 
of time. 

3.5. Dragging gap solitons by a moving lattice 

A problem of considerable interest is a possibil- 
ity of controllable transfer of solitons by a moving 
OL. We have studied this possibility by means of 
direct simulations of Eqs. (fT4)l and pS]) . replacing 
the static OL by V{x) = — Vo cos {2{x -I- vt)), which 
moves at constant velocity —v. Typical results for 
the dragging of GSs taken from the first bandgap are 
displayed in Fig. [51 both for the bosonic model with 
the ordinary cubic nonlinearity, and for the model 
of the BCS superfluid. 

As one might expect, the soliton is dragged in 
a relatively stable fashion at smaller velocities, but 
gets destroyed if v is too large. The quality of the 
transportation may be characterized by the relative 
loss of the number of atoms in the soliton at the end 
of the simulation; these data are collected in Table 
1 (recall that the initial norm of each soliton is 1 , as 
per Eq. P^ ). It is difficult to exactly identify a crit- 
ical value of V at which the moving lattice destroys 
the GS; rather, the transition from the stable drag- 
ging (at w = 2 and 4) to the destruction (at v ~ %) 
is gradual. The results for the bosonic model with 
the quintic nonlinearity (not shown here) are gen- 
erally similar, although in that case the GS appears 
to be more fragile, as its destruction commences at 
velocity v somewhat smaller than in the two models 



Stationary solutions to Eq. 
e" 



(|T8l) are looked for 
as •ip{x,y,t) ~ e^'^*(/)(a;, J/), with real function 4'{x) 
obeying equation 

+ Vo [cos {2x) + cos (2y)] (j) = 0, (37) 

cf. Eq. (pS)) (the small cubic term is dropped here). 
Together with normalization condition pO]) . Eq. 
(j37|) can be derived from the respective Lagrangian, 



L = 



I, 



dxdy <! yu0^ - - {V(j)f 



-G 



(7/3) ,10/3 



2D 



-Vb [cos(2a;) -f cos(2y)] 0^} - m, 



(38) 



cf. Eq. (gT]). To apply the VA, we adopt the 2D Gaus- 
sian ansatz, cf. its ID counterpart (p8)) . as 



(f>ix,y) 



1 
W 



exp 



y 



2W^ 



(39) 



where Af is the 2D norm, that will be actually fixed 
as Af = 1, in accordance with Eq. (PU)) (for the time 
being, M is one of the variational parameters, to- 
gether with fj, and width W, like in the ID setting 
considered above) . An anisotropic ansatz for 2D soli- 
tons, with different widths in the x- and y-directions, 
was considered too. We do not present it here, as nu- 
merical solutions of the variational equations have 
revealed only isotropic solitons. 

The substitution of ansatz (j39p in Lagrangian ( 
yields the effective Lagrangian, cf. Eq. ((29)) ]. 



Left ^fi{M-l)- 



M 



2H/2 



' 257r2/3Ty4/3 



2VoMe- 



-w^ 



(40) 



10 



The first variational equation, dLcs/dj-i = 0, yields, 
as expected, M = 1. Then, equations dLcn/dW = 
dLcff/dM = take the following form: 

-(7/3) 



1 



6g: 



2D 



1 



2M^4 257r2/3 W'lo/a 



1 



/" 



2W^ 



-2V^e 



-W 



3g: 



2^06 

(7/3) 
2D 



-W^ 



1 



(41) 



(42) 



In the linear limit, G2Y) 



57r2/3 M^4/3- 

= 0, Eq. ((4T|) coincides 
with its ID counterpart, Eq. ([32|) . Accordingly, phys- 
ical solutions (W^2 > 0) exist, in the linear limit, for 
^0 > (K)),„in = eVl6 « 0.462, see Eq. ^. Al- 
though the linear Schrodinger equation with a pe- 
riodic potential cannot have localized solutions (in 
any dimension), the result obtained in the linear 
limit makes sense, similar to the situation in the ID 
model: after taking W as the smaller root of Eq. 
((5^ . and then the respective value of ^ from Eq. 
(|42)) . one will obtain a value of fj, at which the fam- 
ily of the GS solutions emerges. In this way, the VA 
predicts the left edge of the first finite bandgap in 
the 2D spectrum (a similar result was obtained in 
|41j . and for the quasi-lD OL potential in the 2D 
setting - also in work [43] , where the VA could accu- 
rately predict the edge of the semi- infinite gap). In 



particular, the VA gives (/i^) 



\ (var) 



-4.26 for 



2D 



Vb = 4, whereas the numerical solution of linearized 



equation ((37)) yields (/Xi)[gfj 



\ (num) 



-4.25 in this 



2D 



case, cf. result (|34|) and its numerical counterpart, 
obtained in the ID case. 

4.2. Numerical results 

Numerical solutions for 2D solitons were obtained 
by running simulations of Eq. (|18p until the solu- 
tions would settle down to stationary real states, 
as done above for ID equation ([T^ : recall this pro- 
cedure produces only stable solutions. The results 
were compared to the predictions of the VA. For the 
purpose of the comparison with BEG, we also in- 
clude findings for the 2D GPE with the usual cubic 
nonlincarity, whose stationary form is 



-G 



(3)j,3 



2D'. 



ti4>+ {1/2) {Cp^^+4>yy) 

+ Vo [cos(2x) 4- cos(2y)] (^ = 0, (43) 

cf. one-dimensional equation (j35p . Solutions to this 
equation are normalized as in Eq. (|20|) . hence the 
respective nonlinearity coefficient is [35] Gj^ = 
2^/2tt [a si a\io) N, where Oho characterizes the tight 
confinement in direction z. 
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Fig. 6. (Color online) Examples of stable 2D GSs, obtained 
in the numerical and variational forms from Eq. I II8I I for 
Vb = 4. The smaller and larger values of G2D correspond to 
the first and second bandgaps, respectively, see Fig. [7] 

Typical examples of the 2D GSs belonging to 
the first and second bandgaps of the 2D spectrum, 
which are displayed in Fig. [6] (and many other ex- 
amples not shown here) demonstrate that the VA 
provides a good fit to the numerically found shapes, 
although, of course, simple Gaussian ansatz (j39p 
does not capture very weak tails attached to the 
central body of the solitons, that become (barely) 
visible in 2D numerical solutions taken near the 
edge of the bandgap. In the bosonic model with the 
cubic nonlinearity, the variational and numerically 
found soliton shapes are very close too (not shown 
here) . 

To quantify the accuracy provided by the VA, we 
calculated the corresponding relative error, 

/ / dxdy I |0„um(x, 2/)p - |0var(a;, y)P | 
err := ■ -, 

J J dxdy\(j)^^^^{x,y)\'^ 

7/3 7/3 

For G2J3 = 2 and Gjq — 10, which arc the cases 

displayed in Fig. [H we have found err = 0.020, 
and 0.079, respectively, which illustrates the accu- 
racy provided by the VA for 2D solitons. With the 
increase of nonlinearity, the shape of the soliton 
slightly deviates from the Gaussian, which gives rise 
to a larger value of error. 

In addition to Eq. p^ with the square-shaped 
OL, we have also considered the same equation with 
the radial-lattice potential (the small cubic term is 
dropped here), 



itp. 



= -ivV.-G(I/^VrV- 



Vo cos (2r) -0 = 0. 

(44) 

Unlike previously studied 2D models with the repul- 
sive cubic nonlinearity and radial potential of the 
Bessel type [33], the potential in Eq. ((331) does not 
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Fig. 7. (Color online) The numerical (continuous lines) and 
variational (strings of symbols) nonlinearity-versus-chemi- 
cal-potential plots for 2D GSs obtained from Eqs. I I18I I and 
1 143 I I. The curve labeled "G = Gip (rad)" additionally dis- 
plays numerical results for radial GSs generated by Eq. 1 144 I I 
with the axisymmetric potential. 



vanish at r -^ oo, hence radial GSs can be pro- 
duced by this equation [32] • The shape of these soh- 
tons (not shown here) is quite similar to that of the 

(7/3) 

GSs supported, at the same values of Gj^ , by the 
square lattice, which suggests that the square OL 
gives rise to nearly isotropic solitons. 

Bearing in mind normalization condition (|20p . 
families of fundamental two-dimensional GSs are 
characterized by dependence Gj^ (p-) (and its 
counterpart, G2Y){fJ.), for bosonic Eq. (|43)) ). similar 
to the ID model. These dependences, as predicted 
by the VA and found from the numerical solutions, 
are displayed in Fig. [7] (the VA for Eq. (|^5)) was 
also based on ansatz (P^ ). Numerical results for 
the family of the radial GSs obtained from Eq. (|44p 
are included too. 

A conclusion is that the VA provides good accu- 
racy in describing the GS families in both models 
with the square-lattice potential, although, as in the 
ID case, the approximation formally predicts that 
the families continue across the narrow Bloch band 
separating the first and second bandgaps, where soli- 
tons cannot exist. Another noteworthy similarity to 
the ID case is that the VA yields a higher accuracy 
for the nonlinearity of power 7/3 than for its cubic 
counterpart. 

The family of 2D GSs was also constructed taking 
into regard the additional cubic term in Eq. ()18p . 
It was found that, as well as in the ID case, this 
term, if taken with physically relevant values of the 
coefficient in front of it (an estimate for which is 
given below), produced a negligible effect. 



4.3. Physical estimates for the two-dimensional 
solitons 

As said above, underlying conditions ([7]) and ([9]) 
definitely hold for the 2D GSs generated by Eq. 
(fT8|) . Proceeding to the estimate for the number of 
atoms in 2D solitons in the BCS supcrfluid, it is 
relevant to mention that no two- (or three-) dimen- 
sional matter-wave solitons of any type, regular or 
gap-mode, have been created, thus far, even in BEG, 
therefore exploring new possibilities for the creation 
of multidimensional matter-wave solitons may be 
quite relevant to the experiment. 

Figure [7] demonstrates that achievable values 
of the normalized nonlinearity strength for 2D 



GSs are G, 



.(7/3; 

^2D 



~ 15. Undoing rescalings 

.X 

(fT7|) and (Uni), we conclude that the corresponding 
number of atoms in the soliton is in the range of 
(A'^2Dsoiiton),„ax ^ 1^^' ^^e rcspcctivc largest quan- 
tum number of the filled states in the transverse 
trapping potential being n ^ 10. This number of 
atoms is sufficient to justify the consideration of 
solitons in terms of the BGS superfluid, and to make 
experimental observation of the solitons possible. 

As for the coefficient in front of the additional 
cubic term in Eq. ((HI), the estimate making use of 
expression ([T9| demonstrates that it may achieve 
values |52d| ~ 0.1, which are higher than \giri\ ~ 
lO"'^ in the ID geometry (see above), but still quite 
small in comparison with G2J3 . 



5. Conclusion 

The objective of this work was to predict the exis- 
tence of quasi-lD and quasi-2D solitons in the BCS 
superfluid formed by a gas of fcrmion atoms, with 
weak attraction between atoms with opposite ori- 
entations of their spins. We considered the experi- 
mentally relevant configuration of the gas trapped 
under the combined action of the ID or 2D opti- 
cal lattice (OL) and tight 2D or ID trap applied in 
the transverse direction(s). The analysis was based 
on the 3D equation for the wave function derived 
from the Lagrangian that included the energy den- 
sity of the weakly coupled BCS superfluid. Then, the 
equation was reduced to effective ID and 2D equa- 
tions, taking into regard the tight transverse con- 
finement. A characteristic feature of the equations is 
the self-repulsive nonlinear term of power 7/3; it was 
also shown that the next-order correction to the en- 
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ergy density of the weakly coupled superfluid gives 
rise to a small additional cubic term in the equa- 
tions. Families of stable one- and two-dimensional 
GSs (gap solitons) were found, by means of the VA 
(variational approximation) and in the numerical 
form, in the first two finite bandgaps of the OL- 
induced spectrum. The VA, based on the Gaussian 
ansatz, provides good accuracy in the description 
of the families of ID and 2D solitons, except for 
very close to edges of the bandgaps, where the GS 
changes its shape from tightly- to a loosely-bound 
one, with undulating tails. In the linear limit, the 
VA accurately predicts borders (i.e., narrow Bloch 
bands) separating different gaps in the spectrum. 
The comparison with the predictions for ID and 2D 
bosonic GSs, supported by the cubic (or quintic) re- 
pulsive nonlinearity, has demonstrated that the VA 
provides a higher accuracy in predicting the solitons 
in the present model with the weaker nonlinearity, of 
power 7/3. Even and odd stable bound states of one- 
dimensional GSs were found too. They may be in- 
terpreted as pairs of fundamental solitons placed in 
adjacent local wells of the lattice. In the 2D case, ad- 
ditionally found solutions are radial GSs supported 
by the radial-lattice potential. 

The quasi-lD and quasi-2D gap solitons, confined 
by the moderately strong transverse potential, may 
contain the number of fermion atoms in the ranges 
of 10"^ — 10^ and 10'^, respectively. The possibility 
of dragging GSs by a moving OL was studied in 
the ID setting. It was demonstrated that there is a 
smooth transition from a regime of stable motion of 
the soliton to its destruction, as the OL's velocity 
increases beyond ^ 1 nun/s. 

Experimental creation of the quasi-lD and quasi- 
2D gap solitons in the BCS superfluid seems quite 
feasible. As concerns further developments of the 
theory, issues of straightforward interest are 3D soli- 
tons, as well as vortex solitons in two dimensions. 
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Appendix: one-dimensional subfundamental soli- 
tons 

The ordinary one-dimensional GPE with the re- 
pulsive cubic nonlinearity and OL potential, which 
corresponds to the first line in Eq. (|35p . gives rise to 
antisymmetric SF (subfundamental) solitons, which 
were found in the second bandgap of the OL- induced 
spectrum in work [45]. They are called so because, 

(3) 

if one uses the notation with fixed Gg^ and arbi- 
trary norm, these solitons have the norm which is 
lower than in the fundamental bosonic GSs, for the 
same chemical potential. A characteristic feature of 
the SF solitons is that two maxima of the density, 
|0(a;)P, are located inside a single cell of the OL po- 
tential (i.e., this antisymmetric soliton as a whole is 
essentially confined to the single cell). The SF soli- 
tons are different from antisymmetric bound states 
of two fundamental solitons, which feature two max- 
ima of I (j){x) I located in different potential wells, see 
Fig. Ha). 

In the ID model derived in this work, i.e., Eq. 
(fT4|) . families of SF solitons can also be found in 
the second bandgap. They are unstable, although 
their instability is weak, as, otherwise, the numerical 
method described above, which is based on the di- 
rect integration of Eq. (fT4|) . would not reveal them. 
Similar to the situation in the GPE, the instabil- 
ity does not completely destroy the SF solitons, but 
rather converts them into fundamental GSs belong- 
ing to the first finite bandgap. 

The SF solitons, as odd solutions to Eq. ([26|) . may 
be represented by the VA based on the accordingly 
modified Gaussian ansatz, 



0(x) 



-1/4 



1^3/2 ^<=^P I 21^2 



(45) 



(its norm is 1 for M = 1). Inserting this ansatz in 
Lagrangian ([?7|) yields 



L = n{M - 1) 



3A'/ 
4W^2 



VaM{l - 2W^)e- 



-W 



25/3r(i3/6) ^(7/3) M'/' 



G 



(5/3)''/''7r5/6 1° 1^2/3 



(46) 



cf. Eq. ([29]) . The first variational equation follow- 
ing from Eq. ([iS]) . dL/dii = 0, gives M = 1, as 
before. The remaining equations, dL/dW = and 
dL/dM = 0, take the following form, cf. Eqs. ([30|) 
and §^: 
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2ii/3r(i3/6) 



9(5/3)^^/^7r5/6 



^,(7/3) ^4/3 



■^ID 



= ^VoW\i 



2W^)e 



-w^ 



(47) 



(7/3) 
ID 






25/3r(13/6) G 

^ (5/3)'^/V5/6 "W^2/3 



(48) 



Note that setting g'\^^^ 



in Eq. (gT]) yields 
an equation for W'^ whicli has two physical roots 
if Vb exceeds a minimum value, (Vo)„,i,-| ~ 0.37, cf. 
expression ([33]) for the minimum value of Vq in the 
VA based on ansatz (1281). The smaller root, if substi- 



tuted in Eq. 

^(7/3) _ 



with G 



(7/3) 
ID 



0, yields (m2)LT^ 



fj, ( GiD '^'' = ) J which, as shown by the dotted curve 

in Fig. [U accurately predicts the lower/left edge of 
the second bandgap. In particular, this approxima- 
tion predicts (Aii)iefT ~ ^■^'^ f™' ^ — 5, while the 
respective value obtained from the numerical solu- 
tion of the Mathieu equation is fi^^'^'{Vo = 5) w 
1.03, cf. variational prediction (|34p for the border 
between the semi-infinite and first finite gaps. 

The comparison of the dependence G\j^ (/i) for 
the SF solitons in the second bandgap, as predicted 
by Eqs. ([47]) and ([48]) . versus its numerically found 
counterpart is presented in Fig. [5Ja). Again, for 
the purpose of comparison, this figure addition- 
ally includes dependences between the nonlinearity 
strength and chemical potential of SF solitons, as 
obtained in the numerical form and by means of 
the VA (that was also based on ansatz (|45|)) for two 
bosonic equations ([35]) . As seen from this figure and 
Fig. [2] the accuracy of the variational prediction 
for the SF soliton families is worse than it was for 
the fundamental GSs. Nevertheless, the prediction 
is still acceptable for the SF solitons in the BCS- 
superfiuid model. In particular, the shape of the 
SF solitons in this model is approximated by the 
VA quite accurately, as shown in Fig. [Sjjb) . On the 
other hand, it is evident from Fig.[8l[a) that the dis- 
crepancy between the VA and numerical results for 
the SF solitons found in the bosonic models is much 
larger, which strongly confirms the above inference 
that the VA works essentially better in the model 
of the BCS superfiuid than in the BEG models, i.e., 
for the nonlinearity of power 7/3 than for the cubic 
and quintic nonlinear terms. 




Fig. 8. (Color online) (a) The same as in Fig. |2] but for ID 
subfundamental solitons, found in the second finite bandgap. 
(b) A typical example of the subfundamental soliton found 
from the numerical solution of Eq. I I26I I. and its variational 
counterpart. 
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